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Abstract 

A  coupled  heat  flow  and  moisture  flow  model  (FROSTB)  was  used  to  simulate 
large  scale  freeze-thaw  experiments  to  assess  its  ability  to  predict  soil  moisture 
conditions  during  freeze  and  thaw.  The  experimental  data  consists  of  temper¬ 
ature  and  soil  moisture  profiles  through  freeze-thaw  cycles  of  a  1  -m  layer  of 
frost-susceptible  silty  sand  over  roughly  2  m  of  gravely  sand.  Two  experimental 
conditions  were  modeled:  1)  where  the  soil  moisture  was  lower  than  specific 
retention  (less  than  1 2%  by  weight)  and  no  water  table  was  present  (dry  case) 
and  2)  where  the  soil  was  fairly  wet  and  the  water  table  was  approximately 
1  m  deep  (wet  case).  During  freezing,  FROSTB  tends  to  predict  ice  contents 
higher  than  those  observed,  which  causes  the  simulated  soil  column  to  thaw 
slower.  During  thawing  the  predicted  moisture  contents  in  the  thawed  soil  were 
close  to  the  measured  values  for  the  wet  case  but  were  always  higher  than  the 
measured  moisture  contents  for  the  dry  case.  Possible  reasons  for  the  dis¬ 
crepancy  are  discussed. 


Cover:  Moisture  distribution  (right)  during  freezing  of  an  unsaturated  soil l 
showing  the  frozen  layer  on  top  underlain  by  a  drier  zone  and  then  the 
saturated  soil  (water  table).  The  moisture  migration  during  freeze-thaw 
is  modeled  by  dividing  it  into  discrete  nodes  and  elements  with  assigned 
material  properties  (left). 

For  conversion  of  SI  metric  units  to  U.S./British  customary  units  of  measurement 
consult  Standard  Practice  for  Use  of  the  International  System  of  Units  (SI),  ASTM 
Standard  E380-89a,  published  by  the  American  Society  for  Testing  and  Mater¬ 
ials,  1916  Race  St.,  Philadelphia,  Pa,  19103. 
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NOMENCLATURE 


Aw,  a 

Gardner  fit  coefficients  for  soil  moisture 

X 

coordinate  (positive  downward) 

ak,  P 

Gardner  fit  coefficients  for  hydraulic 
conductivity 

<l>w 

liquid  water  flow  potential;  three-phase 
system 

E 

phenomenological  calibration  factor  for 
partly  frozen  soil 

0aw 

flow  potential;  air-water-soil  system 

Waw  =  hp) 

S 

gravitational  constant 

^iw 

flow  potential;  ice-water-soil  system 

h 

total  hydraulic  head  ( h  =  hp  +  he) 

ea 

volumetric  air  content 

he 

elevation  head  ( he  =  -x) 

0i 

volumetric  ice  content 

hp 

K 

pressure  head  ( hp  =  u/yw) 

saturated  hydraulic  conductivity  (unfro¬ 

On 

volumetric  unfrozen  water  content  factor 
for  frozen  soil 

zen  soil) 

0o 

porosity 

K¥ 

hydraulic  conductivity  of  partly  frozen 

es 

volumetric  segregated  ice  content 

soil 

0u 

volumetric  water  content  (unfrozen) 

% 

hydraulic  conductivity  (unfrozen  soil) 

Yw 

unit  weight  of  water  (yw  =  gpw) 

t 

time 

Pi 

density  of  ice 

u 

Wgrav 

pore  fluid  pressure 
gravimetric  water  content 

Pw 

density  of  water 
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INTRODUCTION 

This  study  assessed  the  accuracy  of  the  FROSTB 
freeze-thaw  model  for  predicting  soil  moisture 
during  freeze  and  thaw.  The  results  of  simulations 
were  compared  to  a  data  set  compiled  during  a 
vehicle  mobility  project  conducted  in  CRREL's 
Frost  Effects  Research  Facility  (FERF),  a  building  in 
which  full-scale  field  tests  are  conducted  with  soils 
up  to  3.6  m  (12  ft)  deep  (Shoop  et  al.  1991).  FROSTB 
was  used  to  simulate  two  freeze-thaw  cycles:  one 
with  relatively  dry  conditions  and  one  in  which  a 
water  table  was  established  at  0.9-1 .2  m  (3-4  ft) 
below  the  surface .  The  simulation  predictions  were 
compared  with  the  measured  soil  temperature, 
soil  moisture  and  frost  heave.  This  study  empha¬ 
sized  soil  moisture  predictions,  which  were  com¬ 
pared  to  several  data  sets:  moisture  tension,  mois¬ 
ture  contents  from  frozen  core  and  thawed  soil 
samples,  and  water  table  measurements  at  the 
standpipes. 

BACKGROUND  ON  THE  MODEL 

FROSTB  is  a  one-dimensional  coupled  heat  flow 
and  moisture  flow  model  that  computes  frost  heave 
and  thaw  settlement  of  a  pavement  or  soil  profile 
with  time.  It  also  calculates  soil  temperature,  mois¬ 
ture  stress,  water  content,  ice  content  and  density 
through  the  depth  of  the  profile  at  each  time  incre¬ 
ment.  Berg  et  al.  (1980)  originally  developed 
FROSTB  in  a  cooperative  study  funded  by  the 
Corps  of  Engineers,  the  Federal  Highway  Admin¬ 
istration  and  the  Federal  Aviation  Administration. 
Additional  details  beyond  those  described  here  are 
given  in  Guymon  et  al.  (1993).  The  model  assumes 
one-dimensional  vertical  heat  and  moisture  flux 
and  is  based  on  a  numerical  solution  technique 
termed  the  nodal  domain  integration  method.  This 
method  allows  the  same  computer  program  to  be 
used  to  solve  a  problem  by  either  the  finite  element 


Element  Node 

No.  No. 


Actual  Soil  Column  Mathematical 

Representation  Model  Column 


Figure  1.  Example  of  a  soil  profile  divided  into 

nodes  and  finite  elements. 

method,  the  integrated  finite  difference  method  or 
any  other  mass  lumping  numerical  method. 

Figure  1  shows  how  FROSTB  uses  nodes,  which 
are  exact  points,  to  divide  the  column  of  material 
into  horizontal  elements.  Material  properties  are 
assigned  to  the  elements. 

The  program  was  developed  for  solving  prob¬ 
lems  of  seasonal  freezing  and  thawing  of  nonplastic 
soils  and  is  based  on  the  following  primary  assump¬ 
tions,  with  additional  assumptions  reported  in 
Guymon  et  al.  (1993): 

•  Darcy's  law  applies  to  moisture  movement  in 
both  saturated  and  unsaturated  conditions. 

•  The  porous  media  are  nondeformable  as  far  as 
moisture  flux  is  concerned;  i.e.  consolidation  is 
negligible. 


•  All  processes  are  single  valued;  i.e.  hysteresis  is 
not  present  in  relationships  such  as  the  soil- 
water  characteristic  curve. 

•  Water  flux  is  primarily  as  a  liquid;  i.e.  vapor 
flux  is  negligible. 

The  governing  equation  used  in  FROSTB  to  de¬ 
scribe  soil  moisture  flow  is  derived  by  substituting 
the  extended  Darcy  moisture-flow  law  into  the  one¬ 
dimensional  continuity  equation  for  an  incom¬ 
pressible  fluid  flowing  through  porous  media: 


±\K  ^L<*k+-£L 

dx  _  H  dx  dt  pw  dt 


(1) 


where  KH  =  unsaturated  hydraulic  conductivity 
(permeability)  (cm/hr) 
h  =  total  hydraulic  head  (cm  of  water) 
x  =  depth  (cm) 

0U  =  volumetric  unfrozen  water  content  (%) 
Pi  =  density  of  ice  (g/ cm3) 
pw  =  density  of  water  (g/cm3) 

0;  =  volumetric  ice  content  (%) 
t  =  time  (hr). 

The  total  hydraulic  head  h  equals  the  sum  of  the  pore 
pressure  head  ( hp  =  u/yw,  where  u  is  the  pore  water 
pressure  and  yw  is  the  unit  weight  of  water)  plus  the 
elevation  head  ( he  =  -x,  where  x  is  measured  verti¬ 
cally  downward).  The  ice  sink  term,  p;90i/pw3f,  ex¬ 
ists  only  for  a  freezing  or  thawing  zone,  and  in  these 
zones  eq  1  is  coupled  to  the  heat  transport  equation. 
The  ice  sink  term  assumes  that  0;  is  a  continuous 
function  of  time. 

In  FROSTB  the  soil  water  characteristics  are  repre¬ 
sented  using  an  equation  in  the  form  of  Gardner's 
(1958)  function: 

0  = - ^ -  (2) 

u  a  v  7 

Av  frp  +1 


where  0O  =  soil  porosity  (%) 

hp  =  pore  pressure  head  (cm  of  water) 

Aw  =  Gardner's  multiplier  for  the  moisture 
characteristics 

a  =  Gardner's  exponent  for  the  moisture 
characteristics. 

For  each  soil  to  be  modeled,  point  values  of  0U  and  hp 
are  determined  in  a  laboratory  moisture  retention 
test  (Ingersoll  1981).  Equation  2  is  then  fit  to  the  data 
using  a  least-squares  approach  to  determine  the 
best-fit  parameters  Aw  and  a. 

Unsaturated  hydraulic  conductivity  is  also  ap¬ 
proximated  in  FROSTB  using  an  equation  from 
Gardner: 


*»=  \  (3) 

AK\hp\  +1 

where  Kh  =  unsaturated  hydraulic  conductivity 
(cm /hr) 

ks  =  saturated  hydraulic  conductivity  (cm/ 
hr) 

Aj<  =  Gardner's  multiplier  for  hydraulic  con¬ 
ductivity 

p  =  Gardner's  exponent  for  hydraulic  con¬ 
ductivity. 

Point  values  of  and  hp  for  each  soil  are  deter¬ 
mined  in  the  laboratory  by  an  unsaturated  hydrau¬ 
lic  conductivity  test  (Ingersoll  1981),  and  eq  3  is  fit  to 
the  data  using  a  least-squares  approach  to  deter¬ 
mine  the  best-fit  parameters  Ak  and  p. 

Within  the  partially  frozen  zone,  FROSTB  re¬ 
duces  the  unsaturated  hydraulic  conductivity  using 
an  empirical  constant,  termed  the  E-factor,  com¬ 
bined  with  the  ice  content  according  to  the  follow¬ 
ing  equation: 

KF=^P)|-10~£eG  ESi^O  (4) 

where  Kr,  is  the  adjusted  hydraulic  conductivity  in 
a  partially  frozen  element  (cm/hr)  and  E  is  the  em¬ 
pirical  constant  (dimensionless).  The  E-factor  can 
be  set  by  the  user  or  determined  within  the  FROSTB 
program  using  an  empirically  derived  equation 
based  on  the  saturated  hydraulic  conductivity,  ks, 
in  centimeters  per  hour: 

E  =  |(Es-3)2+6.  (5) 

An  example  of  the  function  predictions  relative  to 
calibration  points  is  shown  in  Figure  2.  Additional 
discussion  of  the  E-factor  is  given  in  Guymon  et  al. 
(1993). 

Frost  heave  is  estimated  from  the  total  amount  of 
ice  segregation  in  the  frozen  zone  by 

0s=0i-(0o-0n)  (6) 

where  0S  is  the  volumetric  segregated  ice  content 
(%)  and  0n  is  the  minimum  volumetric  unfrozen 
water  content  (%).  If  0S  is  greater  than  0,  ice  segrega¬ 
tion  has  occurred  and  the  frost  heave  is  computed 
by  multiplying  0S  by  the  element  length.  The  0n 
parameter  establishes  the  pore  water  stress  at  the 
freezing  front  for  the  solution  of  the  moisture  trans¬ 
port  equation.  In  this  study,  0n  was  obtained  by 
assuming  a  moisture  tension  of  -800  cm  of  water 
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and  solving  eq  2.  Thaw  consolidation  from  ice  melt¬ 
ing  is  the  reverse  process  of  that  described  above  for 
ice  segregation.  Upon  thawing,  water  in  excess  of 
the  porosity  is  treated  as  a  source,  forcing  upward 
drainage.  It  is  assumed  that  upon  reaching  the  sur¬ 
face,  the  water  drains  away  laterally. 

To  conduct  the  calculations  described  above, 
FROSTB  requires  the  following  input  for  each  mate¬ 
rial: 

•  Gardner's  coefficients  for  soil  moisture  charac¬ 
teristics; 

•  Gardner's  coefficients  for  hydraulic  conductiv¬ 
ity  characteristics; 

•  Porosity  and  density  of  the  soil; 

•  Thermal  conductivity  and  volumetric  heat  ca¬ 
pacity  of  the  dry  soil;  and 

•  The  E-factor. 

FROSTB  also  requires  the  following  input  for  initial 
and  boundary  conditions: 

•  Element  lengths; 

•  Upper-  and  lower-boundary  pore  water  pres¬ 
sures  with  time; 

•  Upper-  and  lower-boundary  temperatures  with 
time; 

•  Initial  temperature,  pore  pressure  and  ice  con- 


EXPERIMENTAL  VALIDATION 

The  FERF  building  is  temperature  controlled  so 
that  full-scale  freeze-thaw  tests  can  be  conducted 
year-round.  The  building  is  divided  into  test  cells 
measuring  approximately  6.1  m  (20  ft)  square  and 
3.7  m  (12  ft)  deep.  The  mobility  testing  was  con¬ 
ducted  in  six  of  the  northernmost  cells .  Freeze  events 
were  created  by  placing  refrigerated  panels  on  the 
soil  surface.  When  frost  penetrated  the  desired 
amount,  the  panels  were  removed  while  the  ambient 
temperature  in  the  building  was  below  freezing.  The 
building  temperature  remained  below  freezing  while 
a  level  survey  was  conducted  to  determine  the 
amount  of  frost  heave,  and  frozen  cores  were  drilled 
to  retrieve  samples  for  measuring  ice  contents.  Then 
either  the  building  was  heated  or  warm  outside  air 
was  allowed  to  enter  the  building  to  provide  heat  for 
thawing  the  soil  from  the  surface.  During  thaw,  sam¬ 
ples  of  the  thawed  soil  were  extracted  and  tested  for 
density  and  moisture  content. 

The  two  soils  in  the  test  cells,  Lebanon  sand  and 
Pompey  Pit  sand,  have  size  gradations  as  shown  in 
Figure  3.  Both  soils  are  classified  as  sands  with  fines 
(SM)  in  the  Unified  Soil  Classification  System.  Their 


Table  1.  Layer  conditions  in  FERF  mobility  test  cells. 


Depth  range 

Dry  density 

Material 

Condition 

(in.) 

(cm) 

(lb/ft3) 

(Mg/m3) 

Tested  conditions 

Lebanon  sand 

dry 

0-6 

0-15.2 

98.6 

1.551 

wet 

VO 

1 

o 

0-15.2 

101.1 

1.620 

Installed  conditions 

Lebanon  sand 

6-18 

15.2-45.7 

104.8 

1.678 

Lebanon  sand 

18-42 

45.7-106.7 

106.0 

1.697 

Pompey  Pit  sand 

42-138 

106.7-350.5 

123.5 

1.978 
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Figure  4.  Material  depths  and  instrumentation 
locations. 


specific  gravities  are  2.71  for  Lebanon  sand  and  2.74 
for  Pompey  Pit  sand .  They  were  layered  with  1 .07  m 
(42  in.)  of  Lebanon  sand  over  2.44  m  (8  ft)  of  Pompey 
Pit  sand,  which  rests  on  an  underlying  concrete  slab. 
During  placement  the  Pompey  Pit  sand  had  a  rela¬ 
tively  constant  density  and  moisture  content  (Table 
1).  The  Lebanon  sand  had  variable  density  in  three 
layers:  the  bottom  two  layers  resulting  from  the 
placement  preparations,  and  the  top  15.2-cm  (6-in.) 
layer  from  preparations  (watering  or  drying,  tilling 
and  compaction)  for  a  particular  test. 

Three  variables  were  measured  on  a  regular  basis 
to  monitor  the  condition  of  the  test  cells.  Ther¬ 
mistors  measuring  temperatures  were  positioned  at 
2.54-cm  (1-in.)  intervals  down  to  a  depth  of  15.24  cm 
(6  in.),  at  15.24-cm  (6-in.)  intervals  down  to  106.68 
cm  (42  in.)  and  at  30.48-cm  (12-in.)  spacing  down  to 
167.64  cm  (66  in.)  (Fig.  4).  The  thermistors  were 
polled  by  an  automatic  datalogger,  and  readings 
were  recorded  at  2-hour  intervals.  However,  for  the 
wet  case  there  were  a  few  periods  when  data  were 
not  recorded.  Tensiometers,  positioned  as  shown  in 
Figure  4,  were  read  manually  once  per  day  through¬ 
out  the  test  period  to  measure  pore  water  tensions. 


a.  Temperature  (°C).  b.  Tension  (kPa).  The  shaded  area  was  frozen. 

Figure  5.  Measured  data  for  the  dry  freeze  event. 
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a.  Temperature  (°C).  b.  Tension  (kPa).  The  shaded  area  was  frozen. 

Figure  6.  Measured  data  for  the  wet  freeze  event. 


Data  from  frozen  tensiometers  were  disregarded, 
since  they  do  not  work  properly  when  frozen. 
When  a  water  table  was  established  in  the  cells, 
daily  water  levels  were  measured  at  standpipes  at 
the  edge  of  the  test  cells. 

The  two  freeze  events  simulated  with  the  FROSTB 
program  were  characterized  by  a  relatively  dry 
condition  (no  water  table)  and  a  wet  condition 
(with  a  water  table).  The  dry  event  involved  a  cross 
section  that  had  the  installed  density  conditions 
indicated  in  Table  1.  Panels  were  placed  on  the  soil 
surface  from  29  April  to  17  May  1988  (Julian  day 
120-138).  In  the  second  (wet)  event,  the  panels  were 
on  the  section  from  19  September  to  6  November 
1988  (JD  263-311).  The  position  of  the  water  table  in 
the  wet  event  is  a  little  unclear:  the  tensiometers 
indicated  an  average  depth  of  94  cm  (37  in.);  at  the 
same  time,  elevation  measurements  in  the  standpipe 
indicated  a  water  table  depth  of  119  cm  (47  in.).  The 
discrepancy  may  be  because  of  an  uneven  surface 
at  the  test  basin  edge  where  the  standpipes  were 
located  and  variations  in  the  water  table  depth  as 
water  was  added.  The  temperature  and  tension 
data  collected  during  these  dry  and  wet  freeze 
events  are  shown,  respectively,  in  Figures  5  and  6 
(relationships  between  soil  tension  and  water  con¬ 
tent  are  shown  in  Fig.  7). 

Final  heave  was  measured  by  conducting  level 
surveys  at  marked  points  on  the  soil  surface  and 
comparing  them  with  like  measurements  prior  to 
the  freeze  event.  The  moisture  content  of  the  frozen 
cores  was  determined  by  sectioning  the  core  into 
2.54-cm  (1-in.)  pieces  and  determining  gravimetric 
water  content  using  standard  procedures.  The  gravi¬ 
metric  water  content  of  thawed  soil  samples  were 
also  determined.  Fleave  and  water  content  data 
will  be  presented  in  comparison  with  the  predicted 
values. 
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a.  Lebanon  sand. 


b.  Pompey  Pit  sand. 

Figure  7.  Modeled  and  measured  values  of  hydraulic 
conductivity  and  moisture  retention. 

FROSTB  SIMULATIONS— GENERAL  INPUT 
Soil  parameters 

The  various  soil  physical  and  hydraulic  character¬ 
istics  used  to  simulate  material  properties  are  shown 
in  Table  2.  Comparisons  of  the  measured  hydraulic 
data  relative  to  values  predicted  by  the  Gardner's 
equations  are  given  in  Figure  7.  Thermal  properties 
used  in  the  simulations  are  shown  in  Table  3.  The 
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Table  2.  Physical  and  hydraulic  properties  input  to  FROSTB  simulations. 


Lebanon  sand 


Property 

Layer  1 
(0-15  cm) 

Layer  2 
( 15-46  cm) 

Layer  3 
(46-107  cm) 

Pompey 

Pit  sand 

Soil  density  (g/cm3) 

1.551* 

1.678 

1.697 

1.978 

Soil  porosity  (cm3/cm3) 

0.419 

0.419 

0.419 

0.336 

Soil  water  characteristics: 

Aw 

1.962  x  10-5 

1.962  xlO"5 

1.962  x  10"5 

3.7116  x  10“3 

a 

1.975 

1.975 

1.975 

1.268 

Min.  unfrozen  water  cont.  (cm3/cm3) 

0.03602 

0.03602 

0.03602 

0.01787 

Saturated  hydraulic  cond.  (cm/hr) 

1.6 

1.6 

1.6 

5.5 

Permeability  characteristics: 

1.590  x  10-9 

1.590  x  10-9 

1.590  x  10“9 

2.875  x  10-5 

P 

4.623 

4.623 

4.623 

3.806 

*  1.551  g/cm3  for  dry  case,  1.620  g/cm3  for  wet  case. 


value  of  33.5  cal/ cm  hr  °C  was  used  for  soil  solids  in 
all  cases  except  one  (wet  case  6b),  which  used  the 
value  17.0  cal/cm  hr  °C. 


Table  3.  Material  thermal  properties. 


Material 

Specific  heat 
(cal/g  °C) 

Thermal  conductivity 
(cal/cm  hr  °C) 

Water 

1.00 

5.0 

Ice 

0.55 

18.0 

Soil  particles 

0.2 

33.5* 

*  A  value  of  17.0  was  used  for  wet  case  6b. 

Initial  and  boundary  conditions 

The  soil  column  was  simulated  with  the  upper 
boundary  at  the  position  of  the  uppermost  function¬ 
ing  thermistor  and  the  lower  boundary  at  the  bot¬ 
tom  of  the  granular  material,  which  was  located  at 
the  upper  surface  of  a  concrete  slab  beneath  the  test 
cells.  The  total  thickness  of  material  was  350.52  cm 
(11.5  ft).  The  depth  from  the  soil  surface  to  the 
uppermost  thermistor  was  5.08  cm  (2  in.)  for  Cell  10, 
the  dry  case,  and  2.54  cm  (1  in.)  for  Cell  11,  the  wet 
case.  In  all  cases  the  soil  column  was  simulated  with 
99  elements:  2-cm-long  elements  between  the  sur¬ 
face  and  105  cm,  4-cm  elements  between  105  and  225 
cm,  5-cm  elements  between  225  and  245  cm,  and  10- 
cm  elements  from  245  cm  to  the  bottom  of  the 
profile. 

The  upper-boundary  pore  water  pressure  was 
chosen  to  be  computer  generated,  as  follows.  When 
the  profile  is  completely  thawed  and  downward 
vertical  drainage  occurs,  the  surface  pore  water 
boundary  condition  is  modeled  by 


which  means  that  the  velocity  flux  across  this  bound¬ 
ary  is  zero.  The  upper-boundary  condition  is  set  to 
0  cm  of  water  when  the  upper-boundary  tempera¬ 
ture  is  above  0°C  and  frozen  regions  remain  in  the 
column.  When  the  surface  temperature  is  below 


0°C,  a  specified  constant  upper-boundary  pore  pres¬ 
sure  (-300  cm  of  water)  is  used. 

In  general  the  lower-boundary  pore  pressure  con¬ 
dition  of  FROSTB  is  set  by  specifying  discrete  pore 
water  pressures  (tensions)  that  relate  to  the  water 
table  elevation  at  times  when  these  conditions  occur. 
At  intermediate  times,  lower-boundary  pore  water 
pressures  are  linearly  interpolated.  For  all  the  cases 
in  this  study,  we  set  a  constant  lower-boundary  pore 
pressure,  which  produced  a  constant  water-table 
depth  throughout  the  simulation.  Specific  cases  will 
be  discussed  in  a  later  section. 

Input  for  the  upper-boundary  temperature  con¬ 
dition  consists  of  a  set  of  specified  times  and  tem¬ 
peratures  that  are  implemented  as  step  changes .  The 
temperature  input  for  this  study  was  the  measured 
data  at  the  upper  thermistor  for  6-hour  time  incre¬ 
ments.  For  the  wet  case  there  were  periods  with 
missing  data,  and  these  were  estimated  using  a 
linear  interpolation  between  the  closest  available 
data  points. 

When  air  temperatures  are  the  input  values  for 
upper-boundary  temperatures,  FROSTB  adjusts  the 
air  temperature  values  to  represent  the  soil-air  inter¬ 
face  temperatures  using  a  procedure  similar  to  the 
U.S.  Army  Corps  of  Engineers'  n-factor  approach  for 
seasonal  freezing  indices  (Department  of  the  Army 
1966).  An  n-factor  is  defined  as  the  ratio  of  the  sur¬ 
face  index  to  the  air  index  (normally  measured  at  2  m 
above  the  surface),  separately  calculated  for  the  full 
freeze  or  thaw  season.  Because  the  input  values  for 
this  study  were  actual  soil  measurements,  the  n- 
factor  used  was  1.0. 

Bottom-boundary  temperature  conditions  con¬ 
sist  of  a  set  of  times  and  temperatures.  T emperatures 
are  linearly  interpolated  at  intermediate  times.  A 
constant  bottom  temperature  was  specified  in  this 
study  at  a  value  equal  to  the  estimated  value  at  the 
concrete  interface.  The  value  was  calculated  by  ex¬ 
trapolating  the  gradient  between  the  bottom  two 
thermistor  measurements  at  the  start  of  the  simula¬ 
tion. 
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The  initial  conditions  required  to  be  set  for  the 
FROSTB  program  are  the  temperature,  the  pore 
water  pressure  and  the  ice  content  of  each  node. 
The  initial  temperature  condition  was  set  by  inter¬ 
polating  between  the  measured  temperatures  down 
to  the  deepest  measurement  at  167.6  cm  (66  in.)  and 
extrapolating  the  gradient  between  the  bottom  two 
measurements  down  to  the  bottom  of  the  cell.  The 
initial  pore  water  pressure  was  set  by  interpolating 
between  the  measured  tensions  down  to  the  deep¬ 
est  measurement  at  137.2  cm  (54  in.),  beneath  which 
it  was  increased  positively  downwards  with  a  gra¬ 
dient  of  1  cm  of  water  for  each  centimeter  increase 
in  depth.  Simulations  began  before  the  freeze  event 
started,  so  the  initial  volumetric  ice  content  was  set 
at  0.0%  for  all  elements. 

FROSTB  allows  for  a  constant  surcharge  (over¬ 
burden)  pressure  to  simulate  the  pressure  acting 
on  the  top  node  of  the  modeled  column.  This  was 
set  in  the  simulations  to  represent  the  pressure  of 
the  upper  2.54  or  5.08  cm  (1  or  2  in.)  of  soil  that  was 
above  the  shallowest  thermistor. 

The  freezing  point  depression  is  a  constant  value 
that  represents  the  temperature  at  which  water 
freezes.  This  was  set  in  all  cases  to  be  0°C. 

Standard  procedures  used  in  the  FROSTB  calcu¬ 


lations  were  selected  as  follows:  the  fully  implicit 
method  was  used  for  the  moisture  solution,  and  the 
Crank-Nicolson  method  was  used  for  the  heat 
transfer  solution.  Simulations  were  run  with  a  time 
step  of  0.2  hours,  which  is  the  time  at  which  bound¬ 
ary  conditions  are  adjusted;  updates  of  the  thermal 
and  hydraulic  properties  were  set  to  occur  once  per 
hour. 

FROSTB  SIMULATION  RESULTS 
Wet  case 

Four  of  the  simulations  that  were  run  to  model 
the  wet  case  will  be  described.  All  cases  started 
with  the  initial  conditions  shown  in  Figure  8.  Table 
4  summarizes  the  parameters  adjusted  for  these 
simulations  and  the  results. 

In  the  first  case,  referred  to  as  lb,  the  lower¬ 
boundary  pressure  was  set  to  match  the  water  table 
depth  (94  cm)  indicated  by  the  tensiometers,  and 
the  E-factor  (8)  was  calculated  by  the  program.  The 
predicted  results  in  this  case  were  that  the  heave 
would  be  about  twice  that  of  the  measured  range 
and  the  frost  penetration  would  be  slightly  less 
than  measured  (Fig.  9a).  When  the  E-factor  was 
manually  set  to  a  value  of  24,  which  forced  the 


Pressure  (cm  of  water) 


Temperature  (°C) 


-100  0  100  200  300  0  5  10  15 


Figure  8.  Initial  conditions  of  temperature 
and  pore  water  pressure  (wet  case,  Julian 
Day  259)  (10  cm  of  water  =  1  kVa). 


Table  4.  Wet  case  simulation  summary. 

Cell  11/12;  Initial  temps  from  Julian  Day  259;  Upper  boundary  at  2.54-cm  (1-in.)  depth. 


Water  table 


Case 

Lower  pressure 

(cm) 

E-factor 

TKSL* 

Results 

lb 

Match  data  from 

94 

Calc 

33.5 

Heave  higher  than  measured 

tensiometer 

(8) 

Frozen  and  thawed  moisture  high 

4b 

Match  data  from 

94 

Set 

33.5 

Heave/max  frost  match 

tensiometer 

(24) 

Frozen  moisture  high 

Thawed  moisture  closer 

5b 

Match  data  from 

119 

Set 

33.5 

Heave/frost  match 

standpipe 

(24) 

Moisture  high 

Thaw  penetration  quick 

6b 

Match  data  from 

119 

Set 

17.0 

Thaw  penetration  better 

standpipe 

(24) 

Max  frost  too  shallow 

*  Thermal  conductivity  of  soil  particles  (cal/cm  hr  °C). 
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260  280  300  320  340  260  280  300  320  340 

Julian  Day  Julian  Day 


a.  Simulation  lb.  b.  Simulation  4b. 


Julian  Day  Julian  Day 

c.  Simulation  5b.  d.  Simulation  6b. 

Figure  9.  Applied  upper-boundary  temperatures,  predicted  frost  heave  and  frost  penetration  compared  to  measured 
values  for  the  wet  case. 


predicted  heave  to  match  the  measured  (case  4b), 
the  maximum  frost  penetration  also  correlated  well 
(Fig.  9b).  In  a  third  case  (5b)  the  lower-boundary 
pressures  were  changed  to  set  the  water  table  at  the 
level  measured  at  the  standpipe  (119  cm).  This 
resulted  in  slightly  less  heave  and  slightly  more 
frost  penetration  (Fig.  9c).  In  all  of  these  cases  the 
predicted  time  of  complete  thawing  is  shorter  than 
that  measured.  To  match  the  thaw  duration,  an¬ 


other  simulation  was  run  with  a  lower  thermal 
conductivity  of  the  soil  particles  (17.0  cal/ cm  hr°C), 
which  caused  the  end-of-thaw  timing  to  match  but 
reduced  the  predicted  maximum  frost  penetration 
to  less  than  the  measured  amount  (case  6b,  Fig.  9d). 

Contour  plots  of  predicted  temperatures  (Fig. 
10)  show  trends  similar  to  the  frost  penetration 
plots  and  can  be  compared  with  the  measured  data 
shown  in  Figure  6a.  In  general  the  contour  plots  of 
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Figure  10.  Predicted  temperatures  (°C)for  the  wet  case. 


predicted  pore  water  pressure  (Fig.  11)  do  not 
compare  well  with  measured  tension  values  (Fig. 
6b).  The  predicted  values  at  shallow  depths  are 
much  higher  than  those  measured.  The  depth  to  a 
tension  value  of  0  kPa  in  cases  lb  and  4b  do  corre¬ 
late,  since  the  boundary  conditions  were  set  to 
force  this  occurrence.  However,  it  is  encouraging 
that  the  position  of  the  measured  4-kPa  tension  and 
the  predicted  -5-kPa  pressure  contours  are  in  simi¬ 
lar  positions.’1' 

Moisture  data  from  the  frozen  core  and  thawed 
soil  samples  were  acquired  as  gravimetric  water 
contents  (Wgrav).  To  compare  the  gravimetric  data 
with  the  volumetric  values  of  water  (0U)  and  ice  (0;) 
content  predicted  by  FROSTB,  the  predicted  values 
were  converted  to  gravimetric  form  as  follows: 

Wgrav  =  -^  +  09ei  ■  (8) 

dry  density 


*  Note  that  positive  tensions  are  equivalent  to  negative  pore 
pressures. 


Figure  12  compares  the  predicted  values  for  all 
four  wet  cases  with  the  measured  gravimetric  water 
content  from  the  frozen  core  taken  the  day  before 
thaw  began.  The  predicted  values  for  the  frozen  soil 
are  substantially  higher  than  the  measured  values. 

Figure  13  compares  the  predicted  water  content 
with  the  measured  data  from  the  thawed  soil  samples 
for  three  situations  after  thaw  begins.  FROSTB  pre¬ 
dictions  of  water  content  for  the  first  day  of  thaw 
(Fig.  13a)  are  still  much  higher  than  measured.  On 
the  second  day  of  thaw  (Fig.  13b)  the  predictions  are 
much  closer  to  the  measured  data  at  the  shallowest 
depths.  Figure  13c  compares  the  predictions  of  wa¬ 
ter  content  for  the  day  in  which  FROSTB  predicts  a 
thaw  depth  equivalent  to  that  measured  on  the  sec¬ 
ond  day  of  thaw  0ulian  Day  320).  On  that  day 
FROSTB  predicts  water  contents  within  the  thawed 
region  similar  to  those  measured. 

Dry  case 

Three  of  the  simulations  that  were  run  to  model 
the  dry  case  will  be  described.  The  variables  ad¬ 
justed  for  these  simulations,  and  a  short  description 
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Depth  (cm)  Depth  (cm) 


Figure  11.  Predicted  pore  water  pressures  (kPa)for  the  wet  case. 


Figure  12.  Predicted  vs.  measured  gravimetric  water 
content  of  the  frozen  core  taken  on  the  last  day  of  the  wet 
freeze  event  (Julian  Day  313).  The  parameters  for  cases 
lb,  4b,  5b,  6b  are  listed  in  Table  4. 
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a.  First  day  of  thaw. 


b.  Second  day  of  thaw. 


c.  Day  in  which  the  FROSTB-predicted  thaw 
depth  equals  the  depth  on  the  second  day  of  thaw. 


Figure  13.  Predicted  vs.  measured  gravimetric  water  content  of  the  thawed  soil  samples  from  the  wet 
freeze  event. 


of  the  results,  are  given  in  Table  5.  All  three  cases 
started  with  the  initial  temperature  conditions 
shown  in  Figure  14.  Cases  8b  and  9b  started  with 
the  initial  pore  water  pressure  conditions  set  equal 
to  the  measured  conditions  from  the  surface  down 
to  the  position  of  the  deepest  tensiometer  (1.37  m) 
and  then  a  constant  gradient  of  1  cm  of  water  with 
each  centimeter  of  depth  down  to  the  bottom, 
effectively  placing  a  water  table  1.67  m  below  the 
surface,  a  situation  that  was  not  present  in  the 
actual  section  (Fig.  14) .  This  gradient  and  a  constant 


lower-boundary  pore  pressure  were  chosen  so  that 
the  predicted  pressures  would  equal  the  measured 
values  (-30  cm  of  water)  at  the  location  of  the 
deepest  tensiometer.  A  third  case  (10b)  used  initial 
conditions  with  the  measured  pore  water  pres¬ 
sures  down  to  1.37  m  and  a  constant  water  pres¬ 
sure  of  -30  cm  of  water  from  that  point  to  the  bot¬ 
tom,  effectively  moving  the  water  table  to  a  posi¬ 
tion  30  cm  below  the  test  section.  Although  this 
places  the  soil  profile  in  a  state  of  nonequilibrium, 
it  more  accurately  simulates  the  experimental  con- 
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Table  5.  Dry  case  simulation  summary. 

Cell  10;  Initial  temps  from  Julian  Day  110;  Upper  boundary  at  5.08-cm 
(2-in.)  depth. 


Case 

Lower 

pressure 

Water 
table  (cm) 

E- factor 

Results 

8b 

Match  bottom 
tensiometer 
(-30  cm  of  water) 

167 

Calc. 

(8) 

Heave  higher  than  measured 
Moisture  high 

9b 

Match  bottom 

tensiometer 

167 

Set 

(24) 

Heave  and  max  frost  depth  match 
Moisture  high 

10b 

-30  cm  of  water 
at  base 

381 

Set 

(24) 

No  heave 

Frost  slightly  too  deep 

Moisture  still  high 

Pressure  (cm  of  water)  Temperature  (°C) 

-100  0  100  200  2  4  6  8 


Figure  14.  Initial  conditions  of  temperature  and  pore  water  pressure  ( dry 
case,  Julian  Day  110)  (10  cm  of  water  =  1  kPa). 


ditions  at  the  beginning  of  freeze  (i.e.  gravity  drain¬ 
age  of  the  soil  column) .  In  case  10b  a  constant  lower¬ 
boundary  pore  pressure  was  set  to  maintain  the 
water  table  at  30  cm  below  the  test  section  through¬ 
out  the  simulation. 

In  the  first  case  (8b)  the  E-factor  was  set  to  be 
calculated  by  the  program,  which  produced  a  value 
of  8.  As  in  the  wet  case  the  predicted  values  of  heave 
were  higher  than  measured,  and  the  maximum 
predicted  frost  penetration  was  slightly  less  (Fig. 
15a).  The  next  simulation  (9b)  used  a  set  value  of  E 
that  had  been  "calibrated"  earlier  for  the  wet  case. 
Its  predictions  (Fig  15b)  gave  heave  amounts  in  the 
lower  range  of  the  measured  data  and  a  maximum 
frost  penetration  nearly  equal  to  the  measured  depth . 
A  third  simulation  (10b)  used  the  set  E-factor  and 
lower-boundary  pressures  to  simulate  a  water  table 
about  30  cm  below  the  bottom  of  the  cell.  This 
simulation  predicted  no  heave  and  a  maximum 
frost  penetration  deeper  than  the  measured  values 
(Fig.  15c). 

Contour  plots  of  predicted  temperatures  (Fig.  16) 
compare  quite  well  with  the  measured  data  (Fig.  5a). 


Julian  Day 


a.  Simulation  8b. 

Figure  15.  Applied  upper-boundary  temperatures,  pre¬ 
dicted  frost  heave  and  frost  penetration  compared  to 
measured  values  for  the  dry  case. 
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b.  Simulation  9b. 


Julian  Day 


c.  Simulation  10b  (no  predicted  heave  in  this  case). 
Figure  15  (cont'd). 


b.  Simulation  9b. 


Figure  16.  Predicted  temperatures  (°C)for  the  dry  case. 


As  in  the  wet  case  the  predicted  pore  water  pres¬ 
sures  (Fig.  17)  are  much  higher  than  those  mea¬ 
sured  (Fig.  5b). 

Figures  18  and  19  compare  the  moisture  predic¬ 
tions  vs.  the  measured  values  for  the  dry  cases.  The 
predicted  gravimetric  water  contents  for  the  frozen 


core  are  much  higher  than  the  measured  values 
(Fig.  18).  Predictions  for  the  first  day  of  thaw  are 
extremely  high  as  well  (Fig.  19a),  and  even  when 
FROSTB  predicts  a  thaw  depth  equal  to  the  first  day 
of  thaw,  the  predicted  values  are  much  higher  than 
the  measured  values  (Fig.  19b). 
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Depth  (cm)  ^  Depth  (cm)  ^  Depth  (cm) 


a.  Simulation  8b. 


b.  Simulation  9b. 


c.  Simulation  10b. 

Figure  17.  Predicted  pore  water  pressures  (kPa)  for  the  dry 
case. 


Figure  18.  Predicted  vs.  measured  gravimetric  water 
content  of  the  frozen  core  taken  on  the  last  day  of  the 
dry  freeze  event  (Julian  Day  138).  The  parameters 
for  cases  8b,  9b  and  10b  are  listed  in  Table  5. 


a.  First  day  of  thaw. 


b.  Day  in  which  the  FROSTB-predicted  thaw  depth 
equals  the  depth  on  the  first  day  of  thaw. 


Figure  19.  Predicted  vs.  measured  gravimetric  water 
content  of  the  thawed  soil  samples  from  the  dry  freeze 
event. 
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DISCUSSION 

Overall  the  model  predicts  the  frost  penetration 
and  heave  quite  well.  However,  it  tends  to 
overpredict  the  amount  of  ice  formation.  The  addi¬ 
tional  ice  in  the  frozen  soil  then  causes  a  slower 
thaw  than  measured  because  of  the  time  required 
to  thaw  the  excess  ice.  One  reason  for  the  high  pre¬ 
dicted  ice  content  is  the  assumption  that  the  soil 
must  be  100%  saturated  for  frost  heaving  to  occur. 
Moist  soils  will  usually  only  saturate  to  85-95% 
due  to  effective  hydraulic  conductivities  and  por¬ 
osities  (Dirksen  and  Miller  1966).  Additional  air 
entrapment  may  also  occur  during  freezing,  as  air 
is  included  within  ice  crystals  and  ice  lenses,  par¬ 
ticularly  during  rapid  freezing.  Based  on  closed- 
system  freezing  of  unsaturated  soils,  Dirksen  and 
Miller  (1966)  suggested  that  frost  heave  occurs 
when  the  soil  saturation  reaches  90%  rather  than 
100%.  This  is  supported  by  the  water  content  mea¬ 
surements  from  the  wet  case,  which  indicate  an 
average  saturation  of  the  frozen  soil  of  87%  (based 
on  an  average  frozen  density  of 1 .532  g/cm3) .  There¬ 
fore,  modifying  the  model  to  account  for  a  satura¬ 
tion  of  90%  would  more  closely  simulate  the  freez¬ 
ing  process  and  reflect  the  measured  water  con¬ 
tents  and  thaw  progress. 

For  the  dry  case  where  a  water  table  is  not  pres¬ 
ent  yet  measurable  frost  heave  occurs  (average 
heave  =  1.2  cm ),  the  measured  total  water  content 
indicates  an  average  saturation  of  only  40%.  Simu¬ 
lation  9b  generated  frost  heave  within  the  range 
measured,  but  the  model  predicted  an  excess  of  ice 
formation.  Simulation  10b,  the  vertical  drainage 
case  that  more  accurately  simulates  the  initial  soil- 
water  conditions,  had  less  ice  formation  but  no 
frost  heave.  However,  even  simulation  10b  pre¬ 
dicted  that  the  total  water  saturation  of  the  frozen 
soil  would  be  75%,  which  is  still  considerably  higher 
than  the  measured  soil  conditions.  Since  the  pre¬ 
dicted  and  measured  soil  densities  are  nearly  the 
same,  the  reasons  for  overpredicting  ice  content 
and  underpredicting  frost  heave  could  be: 

•  The  laboratory-measured  hydraulic  proper¬ 
ties  of  unfrozen  soil  and  the  field  hydraulic 
properties  of  the  freezing  soil  are  different. 

•  The  physical  process  of  soil  freezing  and  frost 
heave  in  unsaturated  soil  is  not  adequately 
understood. 

•  The  flow  potential  of  the  liquid  water  should 
be  based  on  a  three-phase  system  (air-water- 
ice)  rather  than  a  two-phase  system  (water- 
ice). 

Similar  problems  predicting  soil  moisture  con¬ 


ditions  during  freeze  and  thaw  were  reported  by 
Xia  Xu  et  al.  (1991),  although  they  underpredicted 
rather  than  overpredicted  soil  moisture.  They  sug¬ 
gested  that  the  discrepancy  is  due  to  inaccurate 
measurements  of  soil  hydraulic  properties  and  the 
complications  involved  in  predicting  moisture 
movement  in  a  multiphase  system.  Our  simulation 
uses  laboratory-measured  flow  properties  based 
on  the  soil  density  in  the  experimental  test  bed; 
therefore,  this  is  not  likely  to  be  a  cause  of  error. 
However,  there  are  problems  associated  with  ap¬ 
plying  these  properties  based  on  unsaturated  flow 
(an  air-water-soil  system)  to  a  water-ice-soil  sys¬ 
tem  or  a  water-ice-air-soil  system.  The  model  cur¬ 
rently  accounts  for  changes  in  hydraulic  conduc¬ 
tivity  as  ice  begins  to  form  by  reducing  the  hydrau¬ 
lic  conductivity  values  using  an  E-factor  (eq  4  and 
5).  This  factor  is  empirically  based  and  adjusts  the 
hydraulic  conductivity  for  changes  in  surface  ten¬ 
sions  and  tortuosity  as  ice  forms  in  soil  pores  (Guy- 
mon  et  al.  1993).  Use  of  the  provided  function  to 
calculate  the  E-factor  produced  heave  much  greater 
than  measured  in  this  study.  Further  calibration  of 
the  E-factor  using  frost  susceptibility  test  data 
should  be  examined  to  determine  whether  another 
function  should  be  utilized.  Another  solution  would 
be  to  replace  the  E-factor  with  a  function  that  varies 
pore  water  pressure  with  temperature,  as  has  been 
implemented  in  the  Integrated  Model  (Ly  tton  et  al. 
1990). 

Some  improvement  in  the  soil  moisture  predic¬ 
tions  can  also  be  made  by  adjusting  the  minimum 
unfrozen  water  content.  For  this  study  we  used  a 
value  calculated  by  substituting  a  pore  pressure 
equal  to  -800  cm  of  water  in  eq  2.  Analysis  of  data 
from  an  unfrozen  water  content  test  on  the  Leba¬ 
non  sand  (App.  A)  indicates  that  a  pore  pressure  on 
the  order  of -297  cm  of  water  would  be  developed 
at  a  temperature  of-10°C.  When  a  simulation  was 
run  using  a  minimum  unfrozen  water  content  based 
on  a  pore  pressure  of -297  cm  of  water  (0.167  cm3/ 
cm3),  the  predicted  ice  contents  were  reduced  by 
about  10%  relative  to  those  predicted  with  a  pore 
pressure  of -800  cm  of  water,  but  they  did  not  reach 
the  measured  values.  Changing  this  value  reduces 
the  amount  of  water  being  drawn  from  underlying 
soil  layers  and  produces  less-negative  pore  pres¬ 
sures  in  the  freezing  zone. 

Measurable  frost  heave  in  unsaturated  soils  was 
also  observed  in  other  freeze-thaw  cycles  pro¬ 
duced  in  the  FERF.  Using  the  same  moist  soil  and 
no  water  table,  a  separate  freeze-thaw  cycle  (JD 1 60 
to  201)  resulted  in  an  average  frost  heave  of  1.7  cm 
for  a  frost  depth  of  43  cm  and  total  water  contents 
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in  the  frozen  zone  corresponding  to  an  average 
saturation  of  61 . 1  % .  However,  small  sections  of  the 
frozen  soil  core  indicated  total  water  contents  of  88 
and  89%,  which  would  be  enough  to  generate 
heaving  according  to  our  earlier  arguments.  Based 
on  these  measurements,  it  appears  that  the  heaving 
of  unsaturated  soils  in  a  closed  system  is  isolated 
within  the  soil  profile  and  that  small  lenses  of 
nearly  saturated  soil  and  heave  may  lie  adjacent  to 
relatively  dry  layers  (saturation  of  41-62%).  Some 
of  this  layering  can  be  seen  in  the  profile  of  total 
water  contents  of  the  frozen  core  sampled  at  2 ,5-cm 
increments  (Table  6).  The  water  migration  toward 


Table  6.  Total  water  content  and  saturation  for  the 
freeze  cycle  from  June  8  to  July  19  (JD  160  to  201). 


Water  Water 


Sample 

number 

content 

(%) 

Saturation 

(%)* 

Sample 

number 

content 

(%) 

Saturation 

(%)* 

i 

12.51 

45 

a 

20.15 

73 

2 

11.69 

42 

12 

17.50 

63 

3 

11.84 

43 

13 

17.25 

62 

4 

12.3 

45 

14 

24.54 

89 

5 

12.45 

45 

15 

24.21 

88 

6 

13.74 

50 

16 

18.73 

68 

7 

12.31 

45 

17 

11.28 

41 

8 

16.45 

60 

18 

21.54 

78 

9 

19.65 

71 

19 

17.83 

65 

10 

21.68 

78 

20 

19.48 

71 

*  Based  on  an  average  density  of  frozen  soil  of  1.55  g/cm3. 

Notes:  No  water  table  was  present,  the  frost  heave  was  1.7  cm,  the 
freeze  depth  was  50  cm  and  the  sample  thickness  was  2.5  cm. 


the  freezing  soil  and  the  accompanying  soil  drying 
preceding  the  freezing  front  have  been  carefully 
documented  in  the  laboratory  by  Dirksen  and  Miller 
(1966)  and  Nakano  and  Tice  (1990).  Dirksen  also 
suggested  that  this  process  occurs  on  a  microscopic 
scale  and  therefore  is  not  always  detected  in  our 
macroscale  measurements. 

Most  freeze-thaw  models  assume  that  the  soil 
becomes  saturated  as  it  begins  to  freeze,  form  ice 
lenses  and  heave,  because  this  reduces  the  problem 
to  a  two-phase  system.  Even  for  moist  soil  with  a 
nearby  water  supply,  the  soil  is  not  likely  to  be 
100%  saturated,  as  discussed  earlier,  and  therefore 
even  this  case  is  a  three-phase  system.  Kung  and 
Steenhuis  (1986)  attempted  to  model  heat  and  mois¬ 
ture  transfer  in  an  unsaturated,  partly  frozen  soil, 
but  in  doing  so  they  assumed  that  the  ice  pressure 
is  atmospheric,  a  good  assumption  if  there  is  no 
frost  heave.  When  modeling  frost  heave,  however, 
Miller  (1973)  and  Black  (1991)  supported  the  exist¬ 
ence  and  importance  of  the  pressure  in  the  ice 
phase  based  on  the  geometry  of  the  pore  ice  and 


water  as  predicted  theoretically  by  Miller  (1973)  and 
photographed  by  Colbeck  (1982). 

The  FROSTB  model  currently  calculates  flow  po¬ 
tential  for  the  liquid  water  based  on  the  soil  mois¬ 
ture-tension  curves  measured  in  the  laboratory  and 
assumes  that  this  potential  is  the  same  for  either  an 
air-water  system  (<j>aw)  or  an  ice-water  system  (<j)iw). 
Experiments  by  Koopsman  and  Miller  (1966,  later 
expounded  on  by  Miller  1973  and  Black  and  Tice 
1989)  showed  that  this  is  the  case  for  colloidal  soils 
but  that  qaw  and  <J>iw  are  related  by  the  ratio  of  the 
surface  tensions  (<j)aw  =  2.2  <|>iw)  in  non-colloidal  or 
capillary  soils  such  as  the  granular  Lebanon  sand 
used  in  these  experiments.  Based  on  this  the  liquid 
flow  potential  for  a  three-phase  system  can  be  esti¬ 
mated  by  weighting  the  (j>aw  and  9™  potentials  based 
on  the  percentage  of  liquid  in  contact  with  the  air 
and  water  volumes,  respectively: 


<hv  =(t>a 


0, 


+  <t>r 


vea  +  0iy 


(9) 


where  <]>w  =  liquid  water  flow  potential  for  a  three- 
phase  system 

<j)aw  =  flow  potential  in  an  air-water-soil  sys¬ 
tem 

<j>iw  =  flow  potential  in  an  ice-water-soil  sys¬ 
tem  =  <|>aw/2.2 
0a  =  volumetric  air  content 
0;  =  volumetric  ice  content. 

This  "pseudo"  three-phase  potential  effectively 
reduces  the  liquid  water  flow  potential  in  the  freez¬ 
ing  soil  and  would  therefore  reduce  the  flow  v olume 
and  the  resulting  total  water  contents  in  the  frozen 
zone.  An  added  benefit  of  this  modification  to  the 
flow  potential  is  that  the  controversial  E-factor  may 
have  less  importance. 


CONCLUSIONS 

The  results  of  the  simulations  show  that  FROSTB 
predicts  ice  contents  higher  than  measured  in  the 
FERF  tests.  An  effect  of  the  high  ice  predictions  is  a 
delay  in  thawing,  because  it  takes  longer  to  thaw  the 
excess  ice.  Predictions  of  thawed  moisture  contents 
were  close  to  those  measured  in  the  wet  case,  once 
FROSTB  had  predicted  a  thaw  depth  equal  to  mea¬ 
sured  thaw  depths.  In  the  dry  case  the  predicted 
moisture  contents  were  always  higher  than  mea¬ 
sured,  even  in  case  10b,  where  a  deep  water  table 
was  simulated  and  moisture  predictions  from  equiva¬ 
lent  thaw  depths  were  compared  with  the  measured 
data. 
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Experimental  studies  indicate  that  the  process  of 
freezing  and  frost  heaving  in  unsaturated  soils  (with 
no  water  supply)  occurs  by  heaving  of  nearly  satu¬ 
rated  layers  adjacent  to  relatively  dry  soil  layers  and 
that  this  layering  can  occur  on  a  microscopic  level. 
T otal  water  content  measurements  of  the  frozen  soil 
also  indicate  water  contents  less  than  90%  of  satura¬ 
tion  in  many  cases.  A  possible  solution  to  reducing 
the  amount  of  predicted  ice  would  be  to  allow  ice 
formation  prior  to  saturation  of  an  element,  that  is, 
incorporate  some  air  in  frozen  elements  by  trigger¬ 
ing  frost  heave  at  90%  saturation.  In  addition,  modi¬ 
fying  the  flow  potential  to  account  for  the  presence 
of  three  separate  phases  (eq  9)  effectively  reduces 
the  potential  gradient  for  flow  into  the  frozen  soil. 
The  flow  potential  could  also  be  reduced  by  increas¬ 
ing  the  minimum  volumetric  water  content,  effec¬ 
tively  reducing  the  negative  pore  pressures  in  the 
freezing  zone. 

Another  possible  solution  for  reducing  the  pre¬ 
diction  of  excess  ice  would  be  a  more  rigorous  eval¬ 
uation  of  the  E-factor  function  using  frost-suscepti¬ 
bility  test  data  or  replacement  of  the  E-factor  with  a 
function  that  varies  pore  water  pressure  with  tem¬ 
perature. 
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Unfrozen  Water  Content  (%) 


APPENDIX  A:  UNFROZEN  WATER  CONTENT  DATA  FOR  LEBANON  SAND 


Table  Al.  Warming  data  from  unfrozen 
water  content  test  on  Lebanon  sand.  Four 
samples  tested  with  thawed  gravimetric  wa¬ 
ter  contents  are  shown. 


Temp  Unfrozen  water  content  (%), 

(°C)  gravimetric 


>0 

5.52 

10.87 

15.74 

20.65 

-3.00 

0.32 

0.2 

0.18 

-1.93 

0.63 

0.41 

0.39 

0.37 

-1.46 

0.95 

0.62 

0.59 

0.55 

-0.96 

1.26 

0.82 

0.79 

0.74 

-0.78 

1.26 

1.03 

0.99 

0.74 

-0.59 

1.58 

1.03 

0.99 

0.92 

-0.48 

1.89 

1.24 

1.18 

1.29 

-0.37 

1.89 

1.65 

1.38 

1.29 

-0.29 

2.2 

1.65 

1.58 

1.29 

-0.17 

2.52 

2.06 

1.78 

1.67 

-0.08 

4.09 

2.89 

2.57 

2.04 

-0.04 

4.72 

3.51 

2.96 

2.59 

-0.01 

7.64 

6.32 

5.93 

Table  A2.  Coefficients  for 
modeling  unfrozen  water  con¬ 
tent  of  Lebanon  sand  using 
the  relationship  of  Tice  et  al. 
(1982):  Wu  =  a.  It  |P, where  Wu 
is  gravimetric  unfrozen  wa¬ 
ter  content,  T  is  temperature 
in  °C,  and  a  and  [5  are  coeffi¬ 
cients. 


Sample 

water 

content  Coefficients 


(%) 

a 

P 

5.52 

1.0045 

-0.5629 

10.87 

0.7949 

-0.5069 

15.74 

0.6685 

-0.5297 

20.65 

0.6075 

-0.5235 

Figure  Al.  Warming  data  from  an  unfrozen  water  con¬ 
tent  test  on  Lebanon  sand.  Four  samples  tested  with 
thawed  gravimetric  water  contents  are  shown. 
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1 3.  ABSTRACT  ( Maximum  200  words) 

A  coupled  heat  flow  and  moisture  flow  model  (FROSTB)  was  used  to  simulate  large  scale  freeze-thaw  experi¬ 
ments  to  assess  its  ability  to  predict  soil  moisture  conditions  during  freeze  and  thaw.  The  experimental  data 
consists  of  temperature  and  soil  moisture  profiles  through  freeze-thaw  cycles  of  a  1-m  layer  of  frost-susceptible 
silty  sand  over  roughly  2  m  of  gravely  sand.  Two  experimental  conditions  were  modeled:  1)  where  the  soil 
moisture  was  lower  than  specific  retention  (less  than  12%  by  weight)  and  no  water  table  was  present  (dry  case) 
and  2)  where  the  soil  was  fairly  wet  and  the  water  table  was  approximately  lm  deep  (wet  case).  During  freez¬ 
ing,  FROSTB  tends  to  predict  ice  contents  higher  than  those  observed,  which  causes  the  simulated  soil  column 
to  thaw  slower.  During  thawing  the  predicted  moisture  contents  in  the  thawed  soil  were  close  to  the  measured 
values  for  the  wet  case  but  were  always  higher  than  the  measured  moisture  contents  for  the  dry  case.  Possible 
reasons  for  the  discrepancy  are  discussed. 
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